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ABSTRACT 


Hydrodynamic and semi-kinetic treatments of plasma flow along closed geomagnetic 
fields lines are compared. The hydrodynamic treatment is based on a simplified 16 - 
moment set of transport equations as the equations for the heat flows are not solved; 
the heat flows are treated heuristically. The semi-kinetic treatment is based on a 
particle code. The comparison deals with the distributions of the plasma density, 
flow velocity, and parallel and perpendicular temperatures as obtained from the two 
treatments during the various stages of the flow. In the kinetic treatment, the appro- 
priate boundary condition is the prescription of the velocity distribution functions 
for the particles entering the flux tubes at the ionospheric boundaries; those parti- 
cles leaving the system are determined by the processes occurring in the flux tube. 
The prescribed distributions are half- Maxwellian with temperature T 0 and density 
uq. In the hydrodynamic model, the prescribed boundary conditions are on den- 
sity (no), flow velocity (Vo) and temperature (T 0 ). We found that results from the 
hydrodynamic treatment critically depend on Vo; for early stages of the flow this 
treatment yields results in good agreement with those from the kinetic treatment, 
when Vo = \J A: 7 o/ 27 rm, which is the average velocity of particles moving in a given 
direction for a Maxwellian distribution. During this early stage, the flows developing 
form the conjugate ionospheres show some distinct transitions. For the first hour or 
so, the flows are highly supersonic and penetrate deep into the opposite hemispheres, 
and both hydrodynamics and kinetic treatments yield almost similar features. It 
is found that during this period heatflow effects are negligibly small. When a flow 
penetrates deep into the opposite hemisphere, the kinetic treatment predicts reflec- 
tion and setting up of counterstreaming. In contrast, the hydrodynamic treatment 
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1. Introduction 


In connection with the problem of plasmaspheric refilling, in recent years several 
models for plasma flow along closed magnetic field lines have been developed [ Khaz - 
anov zt al. , 1984; Singh,et aL f 1986; Singh , 1988; Rasmussen and Schunk , 1988; 
Singh , 1991; Guiter and Gombosi , 1990; Wilson, zt al. , 1992]. These models dif- 
fer in complexity in terms of describing the plasma and in including the ionosphere 
as a source of plasma for the refilling. In terms of describing the plasma the most 
contrasting feature of the existing models deals with the hydrodynamic and kinetic 
treatments for the flows, based on plasma fluid equations and a particle-in-cell code, 
respectively. For the purpose of including the ionosphere as a source of plasma for 
the refilling, in most studies the topside ionosphere is replaced by a set of bound- 
ary conditions on the plasma flow, except in Guiter and Gombosi [1990] who have 
included the generation and loss of plasma through chemical reactions in a hydro- 
dynamic model. In this paper, we are mainly concerned with the treatments of the 
plasma with simple sets of boundary conditions on the plasma flow; we compare 
the properties of the plasma flow along closed field lines as given by hydrodynamic 
[Singh, 1992] and semi-kinetic [ Wilson zt al. , 1992] treatments. 

It is generally believed that the kinetic treatment of a plasma is superior to a 
hydrodynamic (fluid) one. The success of a hydrodynamic treatment depends on the 
problem being solved and on the ingenuity of the user in chosing the hierarchy of 
moment equations on which fluid equations are based. In recent years, researchers 
in space physics have used fluid descriptions based on 13-moment [Schunk, 1977; 
Mitchzll and palmadzsso, 1983], and 16-moment [ Barakat and Schunk , 1982; Gan- 
guli and Palmadzsso , 1987; Gombosi and Rasmusszn, 1991; Rorosmezey zt al. , 1992, 
1993] set of transport equations. In developing the moment equations, the ingenuity 
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lies in a series expansion of the plasma distribution function using a biMaxwellian dis- 
tribution function as a base. Therefore, hydrodynamic treatment based on moment 
equations are good as long as the distribution function is close to a biMaxwellian. 
When the distribution function severely departs from a biMaxwellian and involves 
multistreaming of plasma particles, the moment equations are seriously handicapped, 
despite the sophistication incorporated through the use of higher order moment equa- 
tions. 

As mentioned above, recent models for plasmaspheric refilling are based on both 
a kinetic treatment using PIC code and hydrodynamic treatment with varying de- 
gree of sophistication in chosing the hierarchy of the moment equations. In some 
early models only continuity and momentum equations were solved for the ions and 
the electrons were assumed to remain isothermal [Singh et al., 1986; Rasmussen and 
Schunk, 1988; Singh, 1988]. Studies which included temperature equations assumed 
that either the heatflow is given by the collision-dominated thermal conductively 
[Khazanov et al., 1984; Guiter and Gombosi, 1990] or ignored the heatflow com- 
pletely [Singh, 1992; Singh and Chan, 1992]. Neither of these assumptions correctly 
describe the heat transport in the refilling problem [ Singh and Horwitz, 1992]. In 
the collisionless limit of plasma flow during refilling, the usual description of heat- 
flow in terms of Spitzer thermal conductivity breaks down and such a treatment 
overestimates the heatflow. When the heatflow is handled by the 16-moment set of 
equations, the problem remains in the sense that the heatflow equations are valid for 
only relatively small heatflows. When the heatflow becomes large, the validity of such 
equations ceases and numerical instabilities result in computational work. Since, a 
priory it is not known when a large heatflow develops in a model, ad-hoc damping 
mechanisms are included to damp the heatflow, whether it is physically warranted 
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or not [Palmadesso et aL, 1988; Rasmussen, Private Communications]. The above 
problem with the 16-moment transport equations is generic; it arises irrespective of 
the sophistication in numerical techniques employed in solving them [Korosmeze y et 
aL, 1992, 1993]. 

Despite the above difficulties with the hydrodynamic treatment, it has been used 
for practical reasons because it provides simplicity and considerable economy in com- 
putational work and depending on the plasma conditions it can work successfully. 
Therefore, it is advisable to keep in mind the assumptions made in using this treat- 
ment and, if possible it is even better to check the validity of this treatment by 
comparing its prediction against that from a kinetic treatment. Such a comparison 
may reveal when and how a fluid model succeeds. 

The purpose of this paper is to carry out a comparison between models of the 
plasma flows along closed field lines based on kinetic and hydrodynamic treatments. 
The former treatment uses a particle-in-cell (PIC) code for ions [Wilson et aL 1992]. 
The latter one uses transport equations for the flow of mass, momentum, and parallel 
and perpendicular temperatures of ions [Singh, 1992], but the heatflow is treated 
heuristically [Metzler, 1979]. In both the treatments electrons are assumed to obey 
the Boltzmann law. In the present paper, the ionospheric outflows is included by 
imposing a set of boundary conditions on the flow of ions at an altitude of 2000 km. 
The choice of this altitude is primarily due to the existing models [ Wilson et aL, 
1992; Singh, 1992] in which ionospheric loss and generations processes for the plasma 
are not yet included. 

The closed field lines provide the possibility of a variety of flow conditions ranging 
from highly supersonic to subsonic flows as an empty flux tube refills. Futhermore, 
the flows along closed fields lines develop counterstreaming due to interhemispheric 
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plasma flows. Since hydrodynamic treatments are most suspect under counterstream- 
ing situations [ Manheimer et a/., 1976], the comparison carried out here provides a 
useful guide for assessing the validity and usefulness of a hydrodynamic treatment. 

We have found that for the conditions of highly supersonic flows, the two-stream 
hydrodynamic treatment yields flow properties in good agreement with that from 
the semi-kinetic treatment. Demars and Schunk[1991] reported a similar agreement 
based on 16-moment set of equations including heat flows. We find that the bulk 
parameters such as the density, flow velocity and temperatures are in good agreement 
even for the simpler hydrodynamic model when heatflow is included heuristically; the 
reason being simply that when the flow is highly supersonic, the dominant transport 
of heat is through the bulk flow velocity and the transport due to the thermal effects 
is negligibly small. 

When reflection of flows causes counterstreaming, the hydrodynamic treatment 
gives rise to shock formation, which is not seen from the kinetic treatment. How- 
ever, when the counterstreaming flows become subsonic, the hydrodynamic and semi- 
kinetic treatment again produce flow properties in reasonable agreement. 

In order to study plasma flow in space, the plasma treatment must be supple- 
mented by a set of boundary conditions on the flow equation. For the closed field lines, 
the boundary conditions are determined by the top-side ionosphere. The boundary 
conditions along with the demand for plasma at high altitudes produce the flow. The 
ionospheric boundary conditions involves generation and loss of ionospheric plasma 
particle species. Since here our primary goal is in identifying the kinetic and fluid-like 
behaviors of plasma flow and not the supply of plasma from the ionosphere and the 
refilling rate, we have simulated the outflow of ionospheric plasma by imposing a set 
of boundary conditions at an altitude of 2000 km in both the hemispheres. In the 
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semi-kinetic model, the imposed boundary condition is on the velocity distribution 
function of the ions entering the flux tube. It is assumed to be half- Maxwellian. 
The returning particles are self-consistently determined. In the hydrodynamic model 
the boundary conditions are the moments of such a distribution. Since the kinetic 
effect dealing with the returning particles are lost in the hydrodynamic model, the 
hydrodynamic model does not agree with the kinetic model when returning ion flux 
becomes sufficiently large. Can this disagreement be resolved by a more sophisticated 
treatment of the heatflow by using a complete set of 16-moment equations? In order 
to answer this question, further comparative studies are suggested. 

The rest of the paper is organized as follows. The theoretical models are described 
in section 2. The comparison between the results from the two models is carried out 
in section 3. The main conclusions of the paper and their discussion are given in 


section 4. 


2. Theoretical models: 


The semi-kinetic model, which is based on a particle-in-cell code, has been pre- 
viously described for both open [ Wilson, et al., 1990] and closed [Wilson, et al., 
1992] flux tubes. Coulomb collisions are included in the model, the collisions are 
implemented by pairing simulation ions' according to an algorithm which conserves 
energy and momentum [Takizuka and Abe, 197/]. The algorithm yields good ap- 
proximation for the collisions when the collisional relaxation time is shorter than the 
time step in advancing the ion motion. In the hydrodynamic model, we solve the 
plasma transport equations based on 16-moment approximation [e.g., see Ganguii 
and Palmadesso . 1987, and Barakat and Schunk, 1982] 
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where t is time; r is the geocentric distance along the flux tube; s is the distance 
along the tube from the northern ionospheric boundary at A = + ( see Figure 1), 
n, V, Tj| and Tj_ are the number density, flow velocity and parallel and perpendicular 
temperatures of ions in the plasma flow, respectively; 
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< 7 ll and are the heat fluxes along the magnetic field line associated with Xjj and 
Tj_ , respectively; E is the parallel electric field; g\\ is the component of the gravita- 
tional force parallel to the magnetic field, and m and e are the ion mass and charge, 
respectively. The collision terms denoted by [*] c are calculated using Burger s Formu- 
lae [Burger, 1969], which are modified to include flow velocity corrections [Mitchell 
and Palmadesso et al 1983; Ganguli and Palmadesso, 1987] and the correction for 
temperature anisotropy [ Ichimaru et al., 1973; Singh, 1991]. 

We do not solve the heat flow equations, which have proven to be quite trou- 
blesome to solve numerically [ Palmadesso . et al. , 1988]. The difficulty arise for a 
relatively large heatflow, for which the moment equations themselves become invalid. 
Since it is unpredictable in a model when the heatflow may be large, ad-hoc pro- 
cedures are employed to attenuate the heatflow for the numerical stability of the 
models. This has been found to be true irrespective of the numerical techniques 
for solving the equations [ Palmadesso et aL . 1988; Koromezey et al., 1992, 1993; 
Rasmussen, private communication]. 

In this paper instead we have included the effects of heat flow heuristically by 
closely following; the treatments in solar wind studies [e.g. Metzler . et al.. 19(9]. In 
a collisionless plasma the usual picture of heat flow, given by q a — — A a VT a with 
K a as the thermal conductivity, may not be valid because L, the mean free path, 
is » L t = (T-'dT/ds)' 1 , the scale length in the temperature variation. In such a 
collisionless situation, the heat flux can be calculated on physical ground as follows. 
The heat fluxes q ^ and q± across a surface in a single direction in a plasma described 
by a biMaxwellian distribution function with parallel and perpendicular temperatures 
TJl and Tj_, respectively, say along the magnetic field vector, are given by 

qfi = nkT n {kT\\/2irm) l/2 (5) 
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ql = nkT]_{kT\\l2irm) 


( 6 ) 


1/2 

In a uniform plasma for which the distribution function is independent of the parallel 

coordinates, the heatflow at any point is zero because heat flux in a given direction is 

1 

cancelled by the heat flux in the opposite direction. In the presence of a spatial inho- 
mogeniety, the cancellation is not complete, and the heat fluxes <7|| and q± appearing 
in equations (3) and (4) can be heuristically written as [Metzler, et al ., 1979]. 

q a = €T}nkT a V t \\ , (7) 

where the subscript a stands for _L or ||, e = — 1 if dT a /ds > 0 and e — 1 if 
dT a jds < 0. Thus, in the heat flow model adopted here only the sign of the heat flux 
depends on the temperature gradient and not its magnitude. The factor tj determines 
the reduction in the heat flow below the unidirectional fluxes in (5) and (6). Later on 
in this paper we show that 77 in the range say 0.1 - 0.3 yields results in a reasonable 
agreement with the kinetic model, in which heat fluxes appear self-consistently. A 
similar model for heatflow was used by Singh[l 992] for plasma flow along open field 
lines. In both the hydrodynamic and kinetic models adopted here, electric field E is 
calculated by assuming that the electrons obey the Boltzmann law and the condition 
of quasi-neutrality prevails. 

The plasma flow along a closed field line is studied by solving an initial-boundary 
value problem. In the hydrodynamic model, the plasma flows originating from 
the conjugate ionospheres are treated as separate fluids; this treatment is termed 
as a two-stream model [Singh, 1988; Rasmussen and Schunk , 1988; Singh , 1990], 
In both the models, it is assumed that at the initial time (/ = 0) the flux tube 
is highly depleted. The depletion is given by n x — n 0 (s\n A/ sin A 0 ) 6 , with the 
minimum density limited to 10~*no, where no is the density at the ionospheric 
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base A = ±A 0 (Figure 1). Initial flow velocity V(\,t = 0) = 0 and temperatures 
T L {\,t = 0) = 7||(A,f = 0) = T 0 = eV. In the hydrodynamic model, the 
boundary condition for the fluids originating from the northern hemisphere are: 
n n ( A - A 0 ,0 = n*, F n (A = \ 0 ,t) = V 0 , T\\ n {\ 0 , t) = T x {\ 0 ,t) = T 0 ; at the bound- 
ary A = — Xo floating boundary condition are applied. A set of similar boundary 
conditions are used for the fluid originating from the southern hemisphere, but with 
the roles of A = ±A 0 interchanged. In the kinetic model, the boundary conditions on 
ion distribution function /(A, V) are that /(A = A V> 0) and /(A = — A 0 , V < 0) 
are half-maxwellians. with a temperature To. These boundary conditions prescribe 
only the ions entering the flux tubes. The ions leaving the flux tubes are deter- 
mined by the processes occurring inside it. A half-Maxwellian, and not a displaced 
Maxwellian, is chosen because of the following reasons. There is no clear observa- 
tional evidence of supersonic flows along closed field lines at an altitude of 2000 km. 
Futhermore, our calculations show that in about 2 hours the flow in the flux tube 
becomes subsonic nearly all along the flux tube; only for an initial stage of about 2 
hours, supersonic flows are seen. In view of such uncertainties on the the flow velocity 
at 2000 km, a half-Maxwellian serves the purpose of the comparative study. 
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3. Numerical Results 


We compare here the properties of the flow in a flux tube with L — 4 as seen 
from the semi-kinetic and hydrodynamic models. Since the boundary value of the 
flow velocity (V 0 ) and the heat flow factor 77 in the hydrodynamic model are free 
parameters, comparison is performed by varying them over physically reasonable 
ranges. The comparison also deals with the accumulation of plasma in the flux tube 
and the equatorial plasma density. 

3.1 Initial Supersonic Flow 

We recall that the hydrodynamic model is based on two-stream flow in which flows 
originating form the two hemispheres are treated as separate fluids, and the temporal 
evolution of the two streams is separately studied. Likewise, even in the semi-kinetic 
model, the separate identity of the ions originating from the two hemispheres is 
maintained. Figure 2 shows the evolution of the flow originating from the northern 
hemisphere as seen from the semi-kinetic model. This figure gives the phase-space 
density plots in A — Vy plane, where A is the geomagnetic latitude and Vfj is the 
flow velocity along the magnetic field line. The positive and negative values of A 
correspond to the northern and southern hemispheres, respectively. The darkest 
region in the grey-scale plots represents the highest density as indicated by the scale 
on the right-hand side. At time t = .003 hour , the plasma in the tube is essentially the 
initial plasma with a density profile given by n = no [sin A / sin Xq] 6 . At later times this 
plasma expands into the flux tube and it is seen to cross the equator at t = 0.25 hour. 
Along with the expansion, new plasma enters the flux tube at the boundary A = + A 0 . 
It is seen that by the time t = 0.75 hour , the flow has penetrated all the way to the 
opposite boundary at A = — A 0 . It is found that the plasma reaching this boundary 
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is not totally lost, but it is partially reflected back, setting up a counterstreaming 
flow as seen from the plots for t > 1 hour. The reflected flow is seen to reach the 
boundary at A = X 0 by the time t = 2 hours. The plasma flow originating from the 
southern hemisphere shows a similar behavior as shown in figure 2 , with the role of 
boundaries at A = ± A 0 interchanged. It is worth pointing out that after reflections, 
ions merge with the ion stream moving in the opposite direction and they do not 
appear as a separate ion beam. 

The possible consequences of the counterstreaming flow will be discussed later 
on. We now compare the above features of the flow seen from the kinetic model with 
those seen from the hydrodynamic model. Figures 3a to 3d show the comparison for 
t =■ 0.5 hour ; these figures show the distributions of (a) density, (b) flow velocity, (c) 
parallel temperature, and (d) perpendicular temperature. In each panel the curve 
from the kinetic model is labeled, and the curves from the hydrodynamic model 
for three values of the heatflow reduction factor 77 = 0, 77 = 0.05 and 77 = 0.3 are 
indicated by the legend. It is seen from these figures that for most of the flux tube all 
four curves are quite close together, irrespective of the heatflow factor 77 . However, 
near the opposite boundary (A = — A 0 ), the curves from the hydrodynamic model 
tend to diverge from the kinetic model, depending on the value of 77 . This difference 
between the two models is attributable to the fact when the flow begins to slow down 
due to a relative increase in the plasma density, the hydrodynamic model predicts 
an increase in the parallel temperature as clearly seen for A < — 36° in Figure 3c. 
The increase in Tj| enhances the pressure and further slows the flow and enhances 
the density. In the semi-kinetic model, the temperature enhancement does not occur. 
Instead, a counterstreaming develops. This contrast between the two models becomes 
much clearer at later times, for example, at t = 1 hour for which the comparison is 
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shown in Figures 4a to 4d. 

Figure 2 shows that at t = 1 hour the reflected ions set up counterstreaming 
throughout the southern hemisphere (A < 0). Since the hydrodynamic model cannot 
handle the counterstreaming, the reflection process creates a shock, which is clearly 
seen in the density, velocity and temperature plots in Figures 4a to 4c, respectively, 
across the shock indicated by the arrows, density suddenly increases, flow velocity de- 
creases and the parallel temperature also increases. We note that the hydrodynamic 
curves for different values of 77 begin to show some difference among themselves, 
with the curve for 77 = 0.3 being closest to that from the kinetic model. It is worth 
pointing out that the flow velocity in the kinetic model is the average over the coun- 
terstreaming ions. The average velocity is lower than that from the hydrodynamic 
model over the region of the counterstrearmng, but where the counterstreaming has 
not yet occurred (A > 30°) the flow velocity from the two models are generally in 
good agreement. 

The shock first propagates upward to the equator and then downward and reaches 
the ionospheric boundary at A = A 0 at t ~ 2 hours. The propagation of the shock 
in the density profile is shown by the arrows in Figure 5. The transit time of about 
2 hours for the shock is in agreement with the development of the counterstream- 
ing starting in the southern hemisphere and spreading to the northern ionospheric 
boundary by t = 2 hours, (see Figure 2). We find that the heatflow plays only a 
minor role in the motion of the shock; the shock speed is slightly enhanced with 
increased heatflow 77 ; for 77 = 0.3 and 0.05 the shock is already absorbed near the 
boundary A ~ A# , while for 77 = 0, the shock can be still seen in the flux tube at t = 
2 hours . 

After the shock reaches the boundary at A = Ao, the flow in the flux tube becomes 


14 


generally subsonic with respect to the ion- acoustic speed, which is about 10 km/s 
with electron and ion temperatures To = 0.3ev at A = We will discuss the 

subsonic stage after we examine the reason why a shock did not form during the 
early stage ( ~ 1 hour) of the counterstrearmng (Fig. 2) in the semi-kinetic model. 


3.2 Electrostatic shock 

We have just seen that a shock automatically forms in the hydrodynamic model 
as soon as the flow begins to reflect near the opposite boundary. On the other hand, 
the kinetic model does not show the shock formation. Instead, a counterstreaming 
flow develops. We examine this issue in terms of the conditions for shock formation 
and ion velocity distribution function. 

According to the original suggestion of Banks et a/.[19/l], a shock should form 
when supersonic flow from the conjugate hemispheres collide at the equator. The 
flows collide as early as t = 0.25 hour ; the flows from the northern hemisphere can 
be seen from Figure 2 and the corresponding flow from the southern hemisphere is 
the mirror image of this flow with respect to the equator. When the flow begins to 
overlap, the shock should form through the ion-ion instability. The conditions for 
such an instability in a colliding situation are given by 

l -3v tl < V h <MC, and T t > 3T\\ (8) 

where \\ x is the ion thermal velocity, \\ is the ion beam velocity, C s is the ion- 
acoustic speed and M is the mach number, which could be as large as 4 [Forslund 
and Shonk , 1970; Montgomery and Joyce , 1969]. However, it must be mentioned here 
that high critical Mach number M — ^ is determined by the non-linear evolution of 
the electron dynamics including trapping and heating of electrons [e.g., see Singh , 
1988]. For isothermal electrons, as it is assumed in the semi-kinetic model, M = 1 .6 
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[ Tidman and Krai l, 1971]. It is worth pointing out that for large beam velocities, 
oblique ion waves propagating at an angle from the magnetic field are likely to be 
excited. However, the role of such highly oblique waves, which are likely to occur for 
highly supersonic beams, in momentum exchange between interpenetrating beams 
and shock formation is not well understood. 

We examine here the likelihood of the instability occurring from the flow param- 
eters given by the semikinetic model. First we do this exercise for t = 30 minutes 
when the flow has crossed the equator. Figure 6 shows the average flow velocity Vj,, 
the temperature ratio J]j/T e and the ion-acoustic speed C , as function of geomagnetic 
latitude for the flow at t = 30 minutes, shown in Figure 2. Note that the tempera- 
ture ratio is plotted after multiplying it by 10, so that all the plots in Figure 6 can 
utilize the same vertical scale. C , is calculated from C 9 = [k(T c 4- 3T\\)/m] 1 ^ 2 , The 
critical temperature ratio, Tj|/T e = 0.33, for the instability is shown by the segment 
of the thick horizontal line in Figure 6. It is seen that the ions have sufficiently 
cooled down to meet instability condition on the ion temperature over an extended 
equatorial region (|A| < 20°), The flow coming from the opposite hemisphere shows 
a similar feature. The ion- acoustic speed in the equatorial region is about 8 km/s. 
It is seen that over the latitudinal region j A | < 20°, the ion beam velocity is about 
If, ~ 2 C 9 . In the semikinetic and hydrodynamic models discussed here the electrons 
are assumed to obey the Boltzmann law. Therefore, ion beams with such velocities 
are too fast to excite the ion-ion instability and thereby to form shocks in the model. 
Futhermore, it is important to point out that the processes which lead to shock for- 
mation, including the ion-ion instability, are microprocesses, which are suppressed in 
the large-scale models [Singh and Chan , 1993]. If electrons dynamics were rigorously 
included in the model and the associated microprocesses properly resolved, it is likely 
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that the ion-ion interaction would have occurred forming shocks. 

Figure 7 shows the drift velocity 14 , Cs and Tj|/T e at t = 1 hour for the flow 
originating from the northern hemisphere. It is seen that as the ion beam penetrates 
into the opposite hemisphere (A < 0), it gets progressively warmer and the temper- 
ature condition Tj|/X c < 0.3 is not met beyond |A| = 10°. Thus ion instability and 
shock formation are not expected. This indicates that the shock formation in the 
hydrodynamic model (Figures 4 and 5) is an artifact of the model. 

In view of the above discussion in connection with Figures 6 and 7, it emerges 

that on the basis of temperature condition alone, it can be argued that if shocks form, 

they should be during the early stage when the ion beams begin to interpenetrate 

in the equatorial region. During later stages, when the beams penetrate into the 

opposite hemispheres the shock formation is not likely unless some how electrons are 

heated enhancing the temperature ratio T e /T\\. However, as mentioned earlier the 

shock formation in the equatorial region requires a rigorous treatment of electron 

dynamics. In view of the simplified treatment of electrons in the models described 

here, and relatively coarse spatial and temporal resolutions afforded by them the 

issue of equatorial shock formation can not be settled in this paper. 

* 

3.3 Subsonic flow 

After the initial stage of supersonic flows from the conjugate ionospheres, the 
flows become generally subsonic. This is predicted from both the models . Figure 
8 shows the status of the flow at t = 4 hours, from both the hydrodynamic and 
semi-kinetic models. As before, there are three curves from the hydrodynamic model 
which are compared against the curve from the semi-kinetic model. Figure 8b shows 
that the flow velocities obtained for different values of rj from the hydrodynamic 
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model agree with the flow velocity given by the kinetic model. The maximum flow 
velocity of about 5 km/s seen near the boundary A = \o is subsonic with respect to 
the ion-acoustic speed C 5 = 10 km/s. We note that the average flow velocity peaks 
slightly above the boundary in both hydrodynamic and kinetic models. The peaking 
is a consequence of the boundary conditions at A = Ao and the acceleration of ions 
by the pressure and electric field distributions in the close vicinity of the northern 
boundary of the flux tube. 

Figures 8a, 8c, and 8d show that for the subsonic flow the density and temperature 
structures critically depend on the heatflow factor 77. For 77 = 0 . 3 , the hydrodynamic 
model yields results in good agreement with those from the kinetic model. When 
77 becomes too small (77 < 0 . 05 ), the structures in the density and temperature pro- 
files markedly differ from the kinetic model; the density structure shows an extended 
density cavity in the equatorial region, where parallel temperature is relatively high 
[Singh, 1991 ]. Futhermore, for low values of 77 there is density enhancement and cor- 
respondingly a low parallel temperature in the southern hemisphere. When heatflow 
factor is sufficiently large (77 > 0 . 15 ), such structures in n(A) and T||(A) are washed 
away. In a recent paper, Ho et al [ 1993 ] computed the parameter 77 from the semi- 
kinetic model for flows along open flux tubes and it was found to be > 0.1. 

The comparison between the hydrodynamic and kinetic results at t = 12 hours is 
shown in Figures 9a to 9d. The density and temperature structures at this stage are 
qualitatively similar to that at t = 4 hours as shown in Figures 8a to 8d. However it 
is seen that at t = 12 hours , the density and temperature profiles even for 77 = 0.05 
have begun to compare well with that for 77 = 0.3, for which the density distribution 
agrees well with that given by the semikinetic treatment. The discrepancy between 
the densities predicted by the kinetic treatment and the hydrodynamic one for 77 — 
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0.3 is bounded by 15 %, for most part of the flux tube, except near the southern 
boundary A — -A 0 . 

The runs for the comparison between the hydrodynamic and kinetic models were 
carried on until t = J8 hours . For time t > 12 hours, it was found that the hydrody- 
namic model systematically yields densities higher than that given by the semi-kinetic 
model. Figure 10 shows the comparison between the equatorial densities as obtained 
from the two models . This figure shows the temporal evolution of the equatorial den- 
sities found from the kinetic (solid line) and the hydrodynamic (broken line curves) 
models. For the latter model, the densities are plotted for different values of the 
flow velocity ( Vq ) at the boundaries A = ±A <?. We remind ourselves that the results 
from the hydrodynamic model shown in Figures 3, 4, 5, and 8 are for a flow velocity 
V 0 = (kT 0 / 27 rm) 1/2 = 0J9V t . We notice from Figure 10 that for this boundary 
value of the flow velocity, the kinetic and hydrodynamic curves are remarkably close 
for t < 12 hours. This implies that this boundary value of the flow velocity closely 
corresponds to the input flux determined by a half- Maxwellian distribution function, 
which is imposed as boundary condition in the kinetic model. For t > 12 hours, the 
boundary value of V Q = 0 .39V t yields an over-refilling compared to the kinetic model. 
This simply implies that the net influx of ion into the flux tube at the ionospheric 
boundaries steadily decreases in the kinetic model, primarily due to the ions flowing 
out of the flux tube. On the other hand, in the hydrodynamic model, the influx is 
primarily determined by the imposed flow velocity and it remains constant. This is 
demonstrated by comparing the temporal evolution of the total plasma content in 
the flux tube as seen from the two model. 

Figure 11 shows the total content as a function of time. As in Figure 10, for 
the hydrodynamic model the curves are for different values of the imposed velocity 
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at the boundary. It is seen that for Vo = 0.39V t, the hydrodynamic model yields 
nearly the same total plasma content as the kinetic model with nearly the same 
rate of increase in it for t < 12 hours. At later times, the content from the kinetic 
model shows a tendency toward saturation because the rate of increase in the content 
continuously decreases. Even in the hydrodynamic model there is a tendency towards 
the decreasing rate, but the decrease is much slower. This difference in the influx of 
the ions form the two model has a simple explanation. In the kinetic model, some 
of the ions have the liberty to exit the flux tube as they are scattered by Coulomb 
collisions, or as they simply flow out. On the other hand, in the hydrodynamic model 
the plasma entering the flux tube can leave the system only through the opposite 
boundary, where the flow velocity becomes exceedingly small after the shock phase 
(t > 2 hours). This implies that in the hydrodynamic model there is no provision 
for the plasma to leave the system. The slight tendency toward the saturation in 
the hydrodynamic model is due to the changing plasma condition near the boundary 
where the flow originates. As the plasma density near this boundary increases, the 
influx into the flux tends to decrease. 

Figures 10 and 11 also show the equatorial density and the total plasma content 
from the hydrodynamic model for Vo = 0.1 V t and Vo = 0 . It is seen that even for 
V 0 — 0, the equatorial density and the total content are increasing with time and, 
in about 48 hours, they tend to approach the corresponding results from the kinetic 
model. It may sound strange how refilling can occur with a boundary condition of 
V 0 — Q\ The refilling of a flux tube with zero flow velocity as boundary conditions 
in a hydrodynamic model was previously described by Singh et ai , [1986]. 

The ion flux developing near the boundaries from the kinetic model and the 
hydrodynamic model for V 0 = 0.39V^ during the early times [t < 12 hours ) is 1.8 x 10 8 
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ions cm~ 2 s~ l . In the kinetic model, this flux continuously decreases because some 
of the ions entering the flux tube eventually leave. However, in the hydrodynamic 
model the plasma entering the flux tube remains in it, causing the over- refilling for t 
> 12 hours when F 0 = 0.39^. When V 0 = 0 , the hydrodynamic model yields a flux 
of < 10 8 ions cm~ 2 s' 1 and it continuously decreases as the forces (determined by 
density and temperature gradients) on the ions accelerating them into the flux tube 
from the boundary cells diminishes with the refilling. It is worth mentioning that the 
comparison carried out above is based on the simplified boundary conditions in the 
kinetic and fluid treatments and a heuristic treatment of the heatflow. A comparison 
of the plasma treatments without these simplifications will be worthwhile. 

A comparison of plasma distributions in the flux tube at a relatively late time 
(t = 48 hours), as obtained from the two models, is shown in Figures 12 a to 12d. For 
the hydrodynamic model, 77 = 0 . 3 , and the distributions are given for three values of 
the boundary velocity, Vo ~ 0.39V t , 0.lV t , and 0 . Density profiles in Figure 12 a 
show the over-refilling for Vq = 0.39Vt , but when Vo is reduced below 0.lV t , the 
density profiles from the two models disagree near the boundary A = Ap, but away 
from it the agreement considerably improves. The disagreement near the boundary 
is also reflected in the velocity profiles in Figure 12b. Despite the above disagreement 
in the density and velocity profiles near the boundary, the temperature structures 
obtained from the two models are nearly identical. Temperature is nearly isotropic 
( T|| ~ Tx); it rises to about 0.4 eV in the equatorial region from the boundary value 
of 0.3 eV. In the late stage of the refilling, the similarity between the temperature 
profiles, despite the differences in the density and velocity profiles from the two 
models, can be understood by examining the temperature equations (3) and (4) and 
the flow properties. Since the flow velocity is small, the velocity profile has little 
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effect on the temperature profiles. The maximum flow velocity in a localized region 
near the boundary is 2.5 km/s, compared to the thermal velocity of 5.5 km/s and 
ion-acoustic speed of 10 km/s. The density distribution affects the temperature 
distribution through the heatflow terms in equations (3) and (4). In the late stage of 
the flow when the gradients have smoothed out and the densities are relatively large, 
the density distribution also has insignificant effects on the heat flow. 
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4. Conclusion and Discussion 


We have carried out a comparison between semi-kinetic and hydrodynamic mod- 
els for plasma flow along closed magnetic field lines. The comparison has direct 
relevance to the problem of plasmaspheric refilling. It is found that the compari- 
son does not depend only on the plasma physics afforded by the models, but it also 
strongly depends on the boundary condition on the flow velocity. In a kinetic model, 
an appropriate boundary condition is to prescribe the velocity distributions of the 
inflowing ions to be half- Maxwellian for V > 0 at A = X 0 and for V < 0 at A = — \ 0 . 
In the hydrodynamic model, this boundary condition corresponds to the drift veloc- 
ity bo = {kT 0 / S-rnn) 1 ' 2 . A comparison of results from the two models with such 
boundary conditions revealed the following important features of the flows. 

1. When supersonic flows develops in response to a sudden depletion in a 
flux tube, the hydrodynamic and kinetic models yield distribution of den- 
sity, flow velocity and temperatures in generally good agreement. The 
temperature distributions in the region of supersonic flows are found to 
be remarkably similar, showing small effect of the heatflow. It is worth 
pointing out that Demars and Shunk [1991] compared the behavior of 
a highly supersonic plasma flow from a hydrodynamic model based on 
a more complete (16-moment) set of equations with that from a semi- 
kinetic model, demonstrating a good agreement. We have demonstrated 
here that, for a highly supersonic flow, even a much simpler set of hydro- 
dynamic equations are adequate. It is physically explained by the fact 
that the transport of heat in a supersonic flow is dominated by the large 
drift velocity and not by the heatflow process. Mathematically speaking, 
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it implies that the heatflow terms are negligibly small compared to the 
convective terms in the temperature equations. 

2. Both models show reflection of the supersonic flow when it penetrates 
deep into the opposite hemisphere. Since even a two-stream hydrody- 
namic model can not handle the counterstreaming for a given flow, the re- 
flection automatically leads to a shock formation [Rasmussen and Schunk, 
1988; Singh , 1991]. The shock first moves upward toward the equator 
and then downward to the ionospheric boundary. An examination of the 
plasma conditions for shock formation shows that the shock seen in the 
hydrodynamic model is an artifact of the model; the ion beams are found 
to be too warm to excite the ion-ion instability which can subsequently 
produce a shock. The semi-kinetic model shows the development of coun- 
terstreaming for the flow; the counterstreaming advances to the equator 
and downward to ionospheric boundary. It turns out that the transit time 
of the shock all the way to the ionospheric boundary and the time for the 
counterstreaming to spread to this boundary are nearly the same, about 
2 hours. In a previous paper, Singh [1991] reported the shock transit time 
to be about 4 hours, which is in error due to a normalization factor of 2. 
In view of the short transit time of the shock, the shock formation does 
not significantly affect the refilling as evidenced by the comparison of the 
flows from the hydrodynamic and kinetic models for later times. 

Lack of shock formation in the equatorial region, when the ion beams 
begin to interpenetrate [ Banks et a/., 1971] is uncertain in view of the 
spatial and temporal resolutions afforded by a large-scale model and the 
simplicity in handling the electron dynamics by the Boltzmann law. 
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3. After about 2 hours, the flow in each hemisphere becomes subsonic with 
respect to the ion-acoustic speed. This is seen from both the models. 

4. A comparison of the total plasma contents and the equatorial densities 
from the two models indicates a good agreement up to about t ~ 12 hours, 
after which the hydrodynamic model indicates over-refilling of the flux 
tube. The over- refilling is traced to the inability of our hydrodynamic 
model to control the net plasma inflow by the returning particles. The 
inflow is determined by the imposed boundary conditions and the outflow 
of plasma is exceedingly small. On the other hand, in the kinetic model 
the influx gradually decrease due to ions returning from the flux tube, 
showing a tendency toward saturation in the refilling in about 2 days. It 
is worth pointing out that it will be useful to perform a study comparing 
the models based on the kinetic and hydrodynamic treatments by relaxing 
some of the simplifications in terms of boundary conditions and in handling 
of the heat flow in the latter treatment. The boundary conditions can be 
relaxed by including the ionospheric plasma generation processes at low 
altitudes [Gaiter and Gombosi , 1990], 

When the boundary flow velocity in the hydrodynamic model reduces below 
Vq = {kTo / 2irm there is an initial underfilling, but eventually the refilling from 

this model catches up to that given by the semi-kinetic model. For example when 
Vo = 0, the degree of refilling from the two models, in terms of both the equatorial 
density and the total plasma content in the flux tube, becomes approximately the 
same in about 2 days. 

In some previous studies [Singh, et ai, 1986; Rasmussen and Schunk , 1988; Singh, 
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1991], a boundary condition of zero flow velocity was used. It may appear strange 
that a refilling occurs with this boundary condition on the flow velocity. The issue 
is briefly revisited here. 

From the comparison of the plasma contents and the equatorial densities given 
by the models, it is concluded that after about 12 hours, the choice of boundary 
condition in the hydrodynamic model is quite uncertain . In view of this uncertainty, 
the choice of zero-velocity boundary could be useful during the late stage of the 
refilling, it yields under-refilling only near the boundaries, where the density and 
average flow velocity show a discontinuity in the flow. Otherwise, over the rest of the 
flux tube the density and flow velocity are in quite good agreement with those given 
by the kinetic model. 

The hydrodynamic model described here is a two-stream model and includes equa- 
tion for the parallel and perpendicular temperatures. Single-stream hydrodynamic 
models [Singh et al., 1986; Guiter and Gombosi, 1990] suffer from the shortcoming 
that they generate shocks at the equator whether the plasma conditions allow them 
or not. The single- and two-stream models with assumed temperature isotropy suffer 
from the shortcoming that the shock transit time is fairly long and major part of 
the refilling occurs through supersonic flows from the ionospheres [Singh et al.. 1986, 
Rasmussen and Schunk , 1988; Singh, 1991]. This is in contrast to the two-stream 
model with a self-consistent treatment of the temperature anisotropy; this model 
yields evolution from supersonic to subsonic flows at the same time scale as the ki- 
netic model. In this sense, the heuristic treatment of the heatflow described in this 
paper appears to be adequate. This treatment also appears to be adequate even in 
the subsonic stage as long as the flow velocity near the boundary is relatively large 
near the thermal speed, for example, for t < 12 hours in our present calculations. 
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However, when the flow velocity becomes sufficiently low so that a large fraction of 
injected ions in the kinetic model begin to return from the immediate vicinity of the 
boundary, the boundary conditions for the hydrodynamic model diverges from that of 
the kinetic treatment, because this model does not allow for a return flux for a given 
stream. Can this situation be improved by a more rigorous treatment of the heatflow 
and/or by properly including the ionospheric plasma supply [ Guiter and Gombosi , 
1990]? In order to answer this question, it will be useful to compare models based on 
(1) the heuristic heatflow treatment, (2) a more sophisticated treatment of the heat- 
flow using 16-moment set of transport equations, and (3) the semikinetic treatment, 
and all models properly including the analogous ionospheric boundary conditions. If 
the model (1) compares well with the latter ones, the computational effort in using 
transport equations for modeling space plasma will be considerably reduced. 
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Figure Captions: 


Figure 1 . Geometry of a closed magnetic flux tube. The latitudinal angle A and the geocentric 
distance r are shown. The ionospheric boundaries are at s = 0 (A = A 0 ), and s = 

^ max (A — Aq ) 

Figure 2 . Temporal evolution of the flow originating from the northern hemisphere; phase-space 
(s - Vj|) plots are shown. The plot at t = .003 hour nearly shows the initial plasma 
in the flux tube. 

Figure 3. Comparison of flows from semi-kinetic and hydrodynamic models at t = 30 minutes 
For the latter model flows are given for 77 = 0, 0.05 and 0.3: (a) density, (b) flow 
velocity, (c) parallel temperature and (d) perpendicular temperature distributions. 
For most latitudes, the curves are so close together that it is difficult to distinguish 
them. 

Figure 4 . Same as Figure 3, but at t = 1 hour The hydrodynamic curves are distinguished by 
the presence of a shock which is manifested by sudden jumps in density, flow velocity 
and parallel temperature near A ~ —25. Shocks are indicated by the arrows. 

Figure 5 . Propagation of the shock is shown through the jump in the density profiles at (a) 
t =1 hour, (b) t = 1.5 hours and (c) t = 2 hours. For the purpose of comparison 
the curve from the kinetic model and three curves for 77 = 0,0.05, and 0.3 from the 
hydrodynamic model are shown. Shocks are indicated by arrows. 

Figure 6 . Distribution of flow parameters from the kinetic model. The average flow velocity 

V 5 , ion-acoustic speed C s and the temperature ratio Tj|/T e for ions with V|| > 0 are 

%> 

shown for the purpose of instability analysis at t = 30 minutes. The corresponding 
curves for V]| < 0 can be deduced from the symmetry considerations. 
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Figure 7. Same as Figure 6 . but at t = 1 hour. 

Figure 8 . Same as Figure 4. but at t = 4 hours. 

Figure 9. Same as Figure 4. but at t = 12 hours. 

Figure 10 . Temporal evaluation of the equatorial density from the kinetic model and from the 
hydrodynamic model for three values of Vo, 0.39 0.1 V t , and OF 7 -. 

Figure 11 . Temporal evolution of the total plasma content in the flux tube for the kinetic and 
hydrodynamic models. 

Figure 12 . Distribution of (a)density, (b) flow velocity, (c) parallel temperatures and (d) per- 
pendicular temperature at t = 48 hours . All vertical scales are linear in this Figure. 
There are three curves from the hydrodynamic model for Vo = 0, 0.1 \ 1 , and 0.39 V t . 
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